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We demonstrate the possibility of optical nonreciprocal response in a three-mode optomechanical system 
where one mechanical mode is optomechanically coupled to two linearly coupled optical modes simultaneously. 
The optical nonreciprocal behavior is induced by the phase difference between the two optomechanical coupling 
rates which breaks the time-reversal symmetry of the three-mode optomechanical system. Moreover, the three¬ 
mode optomechanical system can also be used as a three-port circulator for two optical and one mechanical 
modes, which we refer to as optomechanical circulator. 


PACS numbers: 42.50.Wk, 42.50.Ex, 07.10.Cm, 11.30.Er 

I. INTRODUCTION 

The fundamental role of nonreciprocal transmission in in¬ 
formation processing has been demonstrated fully by the im¬ 
portant application of electrical diodes in electronic informa¬ 
tion technology with semiconductor p-n junctions. However, 
optical nonreciprocity is constrained by the Lorentz recipro¬ 
cal theorem due to the time-reversal symmetry in linear and 
nonmagnetic media [1]. Traditionally, optical nonreciprocity 
is based on magneto-optical crystals [2] by breaking the time- 
reversal symmetry with the Faraday rotation effect, or opti¬ 
cal nonlinear systems [3] by circumventing the symmetrical 
constraint. Recently, a number of alternative schemes based 
on diverse mechanisms have been proposed, such as spatial- 
symmetry-breaking structures [4], indirect interband photonic 
transitions [5], opto-acoustic effects [6], parity-time symmet¬ 
ric structures [7], and moving systems [8]. Moreover, for the 
potential applications in photonic quantum information pro¬ 
cessing, the realization of nonreciprocal photonic devices with 
the ability to be integrated on a chip and operating on a single¬ 
photon level [9] are desirable features in future. 

With rapidly growing interest as a new class of microscale 
integratable devices, optomechanical systems have shown 
enormous potential for the application in quantum informa¬ 
tion processing [10]. It has already been shown that optome¬ 
chanical systems can be used to induce nonreciprocal effects 
for light [11, 12]. At the beginning, the optical nonreciprocal 
effect is based on the momentum difference between forward 
and backward moving light beams in the optomechanical sys¬ 
tem consisting of an inline Fabry-Perot cavity with one mov¬ 
able mirror and one fixed mirror [11], Subsequently a new 
approach for nonreciprocal optomechanical device was pro¬ 
posed by using strong optomechanical interaction in micror¬ 
ing resonators [12]. The nonreciprocal response is obtained 
for the optomechanical coupling is enhanced in one direction 
and suppressed in the other one by optically pumping the ring 
resonator. In principle, the scheme shown in Ref. [12] can 
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be applied on a single-photon level, in spite of the limitation 
induced by the up-conversion of thermal phonons. 

In this paper, we propose a scheme for optical nonreciproc¬ 
ity in a three-mode optomechanical system, where two optical 
modes are linearly coupled to each other and one mechanical 
mode is optomechanically coupled to the two optical modes 
simultaneously. The two effective optomechanical couplings 
are both enhanced by pumping the two optical modes with 
different external driving fields, respectively. And most cru¬ 
cially, there is a phase difference between the two effective op¬ 
tomechanical couplings, which cannot be absorbed into local 
redefinitions of the operators. Nonreciprocal response of the 
three-mode optomechanical system is induced by this phase 
difference which can be associated with an effective mag¬ 
netic field for the three modes [13, 14]. This mechanism has 
been used in the circuit-QED architecture [13] and phonon 
device [ 14] for breaking time-reversal symmetry, and photon 
or phonon circulator behavior [13, 14] was predicted accord¬ 
ingly. Thus, the present three-mode optomechanical system 
can also be used as a three-port circulator formed by two opti¬ 
cal and one mechanical modes, which we refer to as optome¬ 
chanical circulator. This new type of circulators may serve as 
suitable interfaces for the hybrid network comprised of optical 
(or microwave) and mechanical systems. 

This paper is organized as follows: In Sec. II, the Hamil¬ 
tonian of a three-mode optomechanical system is introduced 
and the spectra of the output fields are obtained formally. The 
optical nonreciprocal response is shown in Sec. Ill and the 
optomechanical circulator behavior is discussed in Sec. IV. Fi¬ 
nally, we draw our conclusions in Sec. V. 

II. MODEL 

We consider a three-mode optomechanical system [15] con¬ 
sisting of two optical modes (a and b, frequencies uj a and 
LOi,) and one mechanical mode (c, frequency u rn ) as shown 
in Fig. 1. The optomechanical coupling rates between the 
optical modes and the mechanical mode are denoted by g a 
and gt,\ two optical modes are linearly coupled mutually at 
rate J and driven by external laser sources with frequencies 
uja.d = 0Jb,d — Wd at rates e a and ej, respectively. The Hamil- 
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FIG. 1: (Color online) Schematic diagram of an optomechanical sys¬ 
tem consisting of two optical modes (a and b ) and one mechanical 
mode (c). The optical modes and the mechanical mode are cou¬ 
pled via radiation pressure, respectively; meanwhile, the two optical 
modes are linearly coupled to each other. 


tonian of the system in the rotating frame at the frequency of 
the driving fields 0 J 4 is 

H = hAad^ a + hA b b^ b + hujmC 1 c + HJ (a^ b + afe^) 
+hg a a l a (c + c^) + Hg b b^b (c + c') 

+ih (e a a^e*^“ + £ b b^e l ^ b — H.c.) , (1) 

where A a = oj a — Ud and A b = ca b — uJd are the detunings 
between the optical modes and the driving fields. Without loss 
of generality, we assume that J, e a and £b are real and </> Q ((f> b ) 
is the phase of laser field coupling to the optical mode a ( b ). 
This kind of Hamiltonian can be realized in the optomechani¬ 
cal system with a membrane in a Fabry-Perot cavity [16], mi¬ 
crotoroid optomechanical cavities [17], optomechanical crys¬ 
tals [18] and electromechanical devices [19]. 

Substituting the Hamiltonian (1) into the Heisenberg equa¬ 
tion and taking into account the damping and corresponding 
noise terms, we get the quantum Langevin equations (QLEs) 
for the operators of the optical and mechanical modes 

-^-a = j-y - i [A a + g a (c + c t )] j a - Ub 

+£ a e l ^ a + v/T^ain, (2) 

= j-y -i [A b + g b (c + J)]}b-iJa 
+£be l<t>b + V^&in, (3) 


= [1 + s„,in (w)] 5 (w + oj'), where the 

term “ 1 ” results from the effect of vacuum noise and ti in (()[,,) 
is the Fourier transform of i> in (t!; n ) for v = a, b, c. 

The mean values of the operators in the steady state can be 
obtained from the nonlinear QLEs (2)-(4) by using factoriza¬ 
tion assumption like ( ca } = (c) (a), and then 

(f + iA') e„e*« - iJ£ b e^ b 
{U) “ (£ + iA' a )(£ + iA 0 + J 2 ’ 

/IA (^ + »A ' a )£be^ b -iJ£ a e^ 

W ” P ~ (^+iA')(?+^)+J 2 ’ 

-i (g a \a\ 2 +g b \/3\ 2 ') 

(C) = ’ 

where A^ = A a +g a (£ + £*) and A' b = A b +g b (£ + £*) are 

the effective detuning including the frequency shifts caused by 
the optomechanical interaction. 

To solve the nonlinear QLEs (2)-(4), we linearize the equa¬ 
tions in the strong driving condition (i.e., e Q j a , £ b 7 b ), 
then the operators are rewritten as the sum of the mean val¬ 
ues and the small quantum fluctuation terms i.e., a = a + 8a, 
b = (3 + 8 b, c = £ + Sc, where 5a -C |a| and 5b -C |/?|. 
Substituting them into the nonlinear QLEs (2)-(4) and keep¬ 
ing only the first-order terms in the small quantum fluctuation 
terms 8 a, 5b, and 8 c, we obtain the linearized QLEs 


(5) 

( 6 ) 

(7) 



~r 5a = 
dt 

(-y - «A' a j 8 a - iG a (5c + 5^) 



iJSb H - yja&ini 

(8) 

dt 

y — iA'^j 5b — iG b (8c + < 5 c^) 



-iJ8a + -v/Tb^in, 

(9) 

■fsc = 

dt 

(-y 1 - 5c - i (G a 8a^ + G*5a ) 



-i (G b 8tf + G* b 5b ) + V 7 ^c in , 

( 10 ) 


where G a = g a u = \G a \e l6a and G b = g b /3 = \G b \e l6b are 

the effective optomechanical coupling rates with phase differ¬ 
ence 9 = 9 b — 9 a . 

For convenience, the linearized QLEs ( 8 )-( 10) can be con¬ 
cisely expressed as 


Jt c = - iUm') c-i (g a a^a + g b tfb)+y/^Ci n . (4) 


dt 


V = -MV + TV in , 


( 11 ) 


Here, 7 a ( 7 b) is the damping rate of the optical mode a 
( b ) and is the mechanical damping rate. ai n , &i n and 
Cj n are the input quantum fields with zero mean values, 
and the spectra of the input quantum fields, s„,i n (w), are 


defined via 


(u/) V m 


U, 


m (w) 5 (ui + oj') and 


where the fluctuation vector V = ( 8 a, 8 b, 5c, 8 a\ 5b* 1 , 5c i ) T , 

the input field vector V in = (a in , b in , c ln , aj n , b\ n , cQ , F = 

diag (y/%, -y/Tm, yfu, y/y>, denotes the damp¬ 

ing matrix and M is the coefficient matrix 
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( 12 ) 


Due to the stability condition, the real parts of all the eigen¬ 
values of matrix M have to be positive. By introducing the 
Fourier transform of the operators 

1 f + °° 

o(uj) = ~^= / o(t)e luJt dt, (13) 

v 27T J — oo 

~ 1 /*+°° 

0 t( W ) = -=/ o* (t) e iult dt, (14) 

V 27T J — oo 

(for any operator o) and using the properties of Fourier trans¬ 
formation, the solution to the linearized QLEs (11) in the fre¬ 
quency domain is 

V (uj) = (M — icoiy 1 rV( n (uj ), (15) 

where I denotes the identity matrix. 

As a consequence of boundary conditions, the relation 
among the input, internal, and output fields is given as the 
following [20] 

?W + Vin = V7 vdv (16) 

for v = a, b, c, and = 7 m . Then the output field vector in 
the frequency domain is 

Vout M = U (uj) V in (uj) , (17) 

where the output field vector V out (uj) is the Fourier transform 

of Vout = ^ a out, ^out, Cout, Oout! ^out) c out) an< ^ 

U (uj) =T (M — iu>I)~ l T — I. (18) 

The spectrum of the output fields is defined by 

s„,out M = J duj' (^vl ut (uj') (w)^ . (19) 

By substituting the expression of VM (uj) [Eq. (17)] into 
Eq. (19), one can obtain [21] 

S'out M = T (lj) Sin M + S v ac (w) . (20) 


Here ^in (Cc?) — (Sa,in (^) 5 &b,in (^) ? Sc,in (cj)) , 

Sout (^) = (Sa,out (^) 5 $b,out (m) 7 &c,out (^)) > 

jn 

Sv&c (uj) = (Sa,vac (uj) , S b ,vac (uj) , Sc,vac (uj)) , and 

/ T aa (uj) T ab (uj) T ac (w) \ 

T(w)= T ba (uj) T b b(uj) T bc (uj) , (21) 

\ T ca (u) T cb (UJ) T cc (uj) J 

where the element (uj) ( i,j = a, 6, c) denotes the scatter¬ 
ing probability that is corresponding to the output field of i 
mode arising from the presence of a single photon (or single 
phonon) in the input field of j mode. The scattering probabil¬ 
ities are given as 


T aa M 

= 

\Uu Ml 2 + 

IM Ml 2 , 

(22) 

Tab (UJ) 

= 

\Ui2 Mi 2 + 

\Ui5 Ml 2 , 

(23) 

T ac (UJ) 

= 

|t/i 3 Ml 2 + 

\Ul6 (uj )| 2 , 

(24) 

T ba M 

= 

\u 2 i Ml 2 + 

\u 2 4 Ml 2 , 

(25) 

T b b (uj) 

= 

\1i22 Ml 2 + 

\u 25 Ml 2 , 

(26) 

T bc M 

= 

\u 2 3 M 2 + 

Me Ml 2 , 

(27) 

Tea M 

= 

IC/31 Ml 2 + 

lt/34 Ml 2 , 

(28) 

T c b M 

= 

M (uj) 2 + 

Ms Ml 2 , 

(29) 

T cc M 

= 

M (uj) 2 + 

Me Ml", 

(30) 

where Ujj (uj) (for i. j = 1, - - - ,6) represents the element 
at the *th row and jth column of the matrix U (uj) given by 


Eq. (18). Su,vac (v = a, b, c ) is the output spectrum contribut¬ 
ing from the input vacuum field, 

Sa,vac (w) = \Uia(uj)\ +|f/i 5 (w)| +|/7ie(w)| ,(31) 

Sfc.vac (w) = \U 2 i(uj)\ + | £^25 (t^) | -\-\U 2 q(uj)\ j (32) 

Sc,vac (tu) = IC /34 (uj )\ 2 + If /35 (w )| 2 + \IJ 3 G (w )| 2 . (33) 


HI. OPTICAL NONRECIPROCITY 


J 


In this and next sections, we numerically evaluate the scat¬ 
tering probabilities to show the possibility of optical nonre¬ 


ciprocal response and optomechanical circulator behavior in 
the three-mode optomechanical system. The optimal parame- 
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FIG. 2: (Color online) Scattering probabilities T ab (to) (black solid lines) and Tf )a (uj) (red dash lines) as functions of the frequency of the 
incoming signal ui for different phase difference: (a) 9 = 0 ; (b) 9 = 7t/ 4; (c) 9 = 7r/2; (d) 6 = 37r/4; (e) 9 = 7r; (f) 9 = 57t/ 4; (g) 9 = 3n/2\ 
(h) 9 = 77 t/4 . The other parameters are A' a = A' b ss = IO7, J — G a = G b e~' e = 'fa /2 =^ b /2 = 7 m /2 = 7/2. 


ters for the observation of optical nonreciprocal response are 
obtained according to the numerical results. The physical ori¬ 
gin for the optical nonreciprocal response and optomechanical 
circulator behavior will be discussed in the next section. 

Scattering probabilities T ab (uj) and T ba (uj) as functions of 
the frequency of the incoming signal uj for different phase 
difference are shown in Fig. 2, where the parameters are 
A; = A' = UJ m = 10 7 , J = G a = G b e~ ie = 7a/2 = 
7 b /2 = 7 m/2 = 7 / 2 . The photon transmission satisfies the 
Lorentz reciprocal theorem [e.g. T ab (uj) = T ba (w)] on the 
condition that 9 = 0 or n. In the regime 0 < 9 < 77 we have 
Tab (uj) < T ba (w); in the regime 7 x < 9 < 2 77 we have 
Tab (uj) > Tba (uj). The optimal optical nonreciprocal re¬ 
sponse is obtained as 9 = 7 t /2 [T a b (uj) ~ 0 and T ba (uj) ~ 1 
at uj = w m ] and 9 = 37 r /2 [ T ab (uj) ~ 1 and T ba (uj) ~ 0 
at ui = uj m \. The condition of G a = G b e~ l9 = 7/2 with 
9a = 9b = g can be obtained approximately by setting the 
amplitudes and the phases of the coupling laser fields as 


7 “m 

~ 0 j 

2 9 

(34) 

7r 

(fra = (jjb — @ ~ ^ ■ 

(35) 


In Fig. 3, the scattering probabilities T ab (uj) and T ba (uj) 
are shown as functions of the frequency of the incoming sig¬ 
nal uj for different effective optomechanical coupling rates G a 
with the parameters: A^ = A' b = uj m = IO7, J = 7 a /2 = 
7 b /2 = 7 m/2 = 7 / 2 , Gb = iG a . It is shown that as the ef¬ 
fective optomechanical coupling is weak ({|G a |, |Gb|} -C 7), 
the scattering probability from b mode to a mode is almost 
the same as the one from a mode to b mode, e.g., T ab (uj) « 
Tb a (uj). With the enhancement of the effective optomechani¬ 
cal coupling rates, the optical nonreciprocal response becomes 
obvious and gets to the optimal effect at about G a = O.57. 

In Fig. 4, we plot the scattering probabilities T ab (uj) and 
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FIG. 3 : (Color online) Scattering probabilities T ab (uj) (black solid 
line) and T ba (uj) (red dash line) as functions of the frequency of the 
incoming signal a i for different effective optomechanical coupling 
rates: (a) G a = O.O57; (b) G a = O.257; (c) G a = O.57; (d) G a = 
7. The other parameters are A' a = A' b = uj m = IO7, J = 7a/2 = 
7b/2 = 7 m /2 = 7/2. Here we fix Gb = iG a which corresponds to 
the case of 9 = n/2. 


Tba (uj) for different mechanical damping rates with the 
parameters: A' a = A' h = uj m = IO 7 , J = G a = —iG b = 
7 a/2 = 7 b /2 = 7 / 2 . It is shown that as the mechanical damp¬ 
ing rates is much smaller than the optical damping rate 
(7m 7 ), the photon scattering probabilities are almost the 

same for the two directions, e.g. T ab (uj) ~ T ba (uj). With the 
increase of the mechanical damping rate, the optical nonre¬ 
ciprocal response becomes obvious and achieves the optimal 
effect for y m « 7 . 
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FIG. 4 : (Color online) Scattering probabilities T a b (cu) (black solid 
lines) and T ba (uj) (red dash lines) as functions of the frequency of 
the incoming signal co for different mechanical damping rates: (a) 
7m = O.OI7; (b) 7m = 0.27; (c) 7m = 71 (d) 7 m = 27. The 
other parameters are A' a = A' b = uj m = IO7, J = G a = — iGb = 
7 a /2 = 76/2 = 7/2. 


IV. OPTOMECHANICAL CIRCULATOR 


As done in most of the studies on optomechanical systems, 
the signal input and/or output from the mechanical mode is not 
considered in last section. With the development of phonon- 
based system, phonon is another useful media for quantum in¬ 
formation processing [14]. In this section, we assume that the 
mechanical mode is coupled to a continuous mode of phonon 
waveguide and the phonons can be input and output through 
the phonon waveguide. The scattering of both photons and 
phonons in the three-mode optomechanical system is consid¬ 
ered in the following. 

Using Eqs. (22)-(30), we now show the numerical results 
of all the scattering probabilities (nine elements) in Eqs. (21). 
As shown in Fig. 5, the three-mode optomechanical system 
shows circulator behavior: when 6 = ir/2, we have T ba (w) ~ 
T c b (w) ~ T ac (oj ) ft: 1 and the other scattering probabilities 
equal to zero at w = oj m as shown in Figs. 5 (a), (c) and (e); 
when 6 = 37t/2, we have T ca ( ui ) ~ T ab (w) ~ T bc (w) ~ 1 
and the other scattering probabilities equal to zero at w = uj rn 
as shown in Figs. 5 (b), (d) and (f). That is to say the signal 
is transferred from one mode to another either counterclock- 
wisely (a —> b —> c —> a) or clockwisely (a —> c —> b —> a), 
depending on the relative phase 9 = 7t/ 2 or 37 t/ 2 as shown in 
Fig. 6. 

The scattering matrix for the optomechanical circulator 
in three-mode optomechanical systems can be obtained an¬ 
alytically, similar to the case for photon and phonon circu¬ 
lators in Refs. [13, 14]. We assume that ~ A 
{J, |G a | , |G&| ,7a,7fe>7m}, then Eqs. (8)-(10) can be simpli- 



FIG. 5 : (Color online) Scattering probabilities T va (ui) [(a) and 
(b)], T vb {to) [(c) and (d)] and T vc (w) [(e) and (f)] (v = a, b, c ) 
as functions of the frequency of the incoming signal u> for differ¬ 
ent phase difference: (a), (c) and (e) 9 = 7r/2; (b), (d) and (f) 
8 = 3n/2. The other parameters are A' a = A' b = u> m = IO7, 
J — G a — G b e~ ie = 7o /2 = 76/2 = 7m /2 = 7/2. 



FIG. 6: (Color online) Schematic diagram of a three-mode optome¬ 
chanical circulator. 


fled by rotating wave approximation as 
^ 5a = (—y - *A' a ) Sa-iGaSc-iJSb+y/^ai n , (36) 

y Sb = (-y - iA'hJ 5b - iG b 5c - iJSa + (37) 

y<5c = ^-^y - iiOrr^j 5c-iG* a 8a-iG* b 5b+y/ 7wc in . (38) 
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Thus the effective Hamiltonian after linearization takes the 
form (effective system is shown in Fig. 6 ) 

H e g = HA^a 1 a + HA' b b^ b + hu> m c^ c (39) 

+h (Jtfa + \G a \e i0a a^c + \G b \e~ ieb Jb + H.c.) . 

The necessary and sufficient condition of time-reversal sym¬ 
metry is the gauge-invariant phase sum is an integral multiple 
of 7 T [13]. That is, 


9 b — 9 a = 9 = mr 


(40) 


for real J, where n is an integral number [also see the nu¬ 
merical results as shown in Figs. 2(a) and 2(e)]. The opti¬ 
cal nonreciprocal response and the optomechanical circulator 
behavior are induced by breaking the time-reversal symme¬ 
try (i.e., 9 b — 9 a ^ mr), and the optimal effect is realized at 
the halfway between the time-reversal symmetric points (i.e., 
9 a = 0, 9 b = 9 = 7 r /2 or 3tt/2) as shown in Fig. 5. 

Now we will derive the scattering matrix of the optome¬ 
chanical circulator behavior analytically from the simplified 
linearized QLEs [Eqs. (36)-(38)]. Let us transform the lin¬ 
earized equations into the frequency domain, 


where 


M' = 



( \ 

( 

' — iujl) 

" = 1 

-y/Tb&in 


[Sc) 

V v/TmGn 

l+*A' a 

iJ 

iG a 

iJ 


lG b 

iG* a 

iGt 



(41) 


(42) 


In the conditions for the optimal optomechanical circulator, 

i.e. w = A' a = A' b = u m and J = G a = G b e~ l9 = 7 Q /2 = 
7 b / 2 = 7 m /2 = 7 / 2 , we have 


7 

2 





( 6 ~\ 

6 b 

\6c) 


V70in \ 

VTfcn 

\/7 c in / 


(43) 


By choosing 9 = 7 r/ 2 , the scattering matrix is given through 




(44) 


= 37 t/ 2 , we can get the scattering matrix 


(45) 



Equation (44) shows clearly a perfect circulator with the sig¬ 
nal transferring counterclockwisely (a —> b —> c —> a) for 
9 = 7t/ 2 and Eq. (45) also describes an ideal circulator but 
with the signal transferring clockwisely (a —> c —> b —> a) for 


x to 3 x 10 3 



FIG. 7: (Color online) The vacuum noise spectrum s v , v ac (tu) (v = 
a, b, c) as a function of the frequency of the incoming signal lj for 
different phase difference: (a) 9 = 7 r/ 2 ; (b) 9 = 3iv/2. The other 
parameters are A' a = A' b = uom = IO 7 , J = G a = G(,e“* 6 = 
7 <j /2 = 7&/2 = 7 m /2 = 7 / 2 . 


9 = 37t/ 2. These agree well with the numerical results shown 
in Fig. 5. 

Finally, we discuss the effects of the vacuum noise spec¬ 
trum s„ jVac (w) given by Eqs. (31)-(33). The vacuum noise 
spectrum s„ jVac (w) ( v = a, b, c) as a function of the frequency 
of the incoming signal ui is shown in Fig. 7. The effects of the 
vacuum noises are so small that they are insignificant even for 
the input signals of single-photon (single-phonon) level (about 
0.2% at w = uj m = IO 7 ). The physical origin of the vacuum 
noise in the output spectrum is the anti-rotating-wave interac¬ 
tions between the optical and the mechanical modes [included 
in Eq. (11)]. The suppression of the vacuum noise for b mode 
(a mode) at ui = as 9 = 7 r /2 (9 = 3n/2) is the conse¬ 
quence of the rotating-wave approximation for the interaction 
between the two optical modes. 


V. CONCLUSIONS 

In summary, we have shown the optical nonreciprocity in 
a three-mode optomechanical system. We demonstrated that 
the nonreciprocal response is enabled by tuning the phase dif¬ 
ference between the optomechanical coupling rates to induce 
the time-reversal symmetry breaking of the system. Then we 
show that the three-mode optomechanical system can also be 
used as a three-port optomechanical circulator for two optical 
modes and one mechanical mode. Further, we note that the 
three-mode optomechanical system can work in the single¬ 
photon level and be integrated into a chip. The three-port op¬ 
tomechanical circulator might eventually provide the basis for 
applications on quantum information processing or quantum 
simulation [ 22 ]. 


Note added 

In the preparation of this work, we became aware of a re¬ 
lated paper by Metelmann and Clerk [23]. 
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